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In this paper we give a theoretical description of resonant coupling between two collective exci- 
tations of a Bose condensed gas (BEC) on, or close, to a second harmonic resonance. Using analytic 
expressions for the quasi-particle wavefunctions we show that the coupling between quadrupole 
modes is strong, leading to a coupling time of a few milliseconds (for a TOP trap with radial fre- 
quency ~ 100 Hz and ~ 10 4 atoms). Using the hydrodynamic approximation, we derive analytic 
expression for the coupling matrix element. These can be used with an effective Hamiltonian (that 
we also derive) to describe the dynamics of the coupling process and the associated squeezing effects. 

(N| \ PACS numbers: 03.75.Fi, 05.45.-a, 42.65.Ky 
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I. INTRODUCTION 



In two recent experiments we observed resonant coupling between the low-energy modes of oscillation in a 
Bose-condensed gas. In the first experiment Q we excited an even parity quadrupole mode (the m = low-lying 
mode) and observed transfer of energy to a mode at twice the original frequency (the m = high- lying mode). 
The oscillations at the second harmonic were observed as soon as the excitation period ended and stayed constant 
in amplitude. This indicates strong coupling between the modes so that energy is transferred between them at a 
rate comparable to the mode oscillation frequency of a few hundred Hz, i.e. this an allowed transition between the 
vibrational modes. In contrast, the coupling between a scissors mode and a mode at half the initial frequency was 
7— I , found to be a much slower process B. This paper shows that the simple downconversion process is forbidden, i.e. 

the matrix element for the direct conversion of one quantum of the higher-lying scissors mode into two quanta of a 
lower-lying mode, is zero. This means that some more complicated process is required to explain the experimental 
£^ ' results. We also show how to calculate the coupling rates between various modes analytically. For resonant coupling 
ly-j , between the quadrupole excitations we present a simple expression for the radial integrand of the matrix element 
f^) 1 that shows that the coupling mostly takes place in the boundary regions of the condensate. Finally, we show that the 
coupling is well described by a simple Hamiltonian which can be used for quantitative studies of the squeezing effects 
1 related to the harmonic generation processes. 

The paper is structured as follows: Section [n] presents the nonlinear Schrodinger equation (NLSE) and the derivation 
of the Bogoliubov-de Gennes (BdG) equation from the many-body Hamiltonian. These equations form the basis of the 



following sections. In section HI we summarize the derivation of solutions to the BdG equations in the hydrodynamic 
limit following the method given in The assumptio ns a nd approximations made in that derivation are important 
for understanding the calculations in section |v|. Section IV gives the derivation of the Hamiltonian describing second- 
^ ■ harmonic generation SHG or degenerate down-conversion from the NLSE, closely following the approach given in [0. 
The coupling matrix elements governing the nonlinear processes are calculated in section 0. A simple expression is 
found for resonant coupling and the results are compared to an exact numerical calculation. We show that symmetry 
arguments forbid the direct down-conversion of the scissors mode and discuss our results with respect to two recent 
jJJ 1 experiments Jl]]^] by our group. 
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II. CONDENSATE EXCITATIONS 



Our treatment of the coupling between the modes starts with the Gross-Pitaevskii equation (GPE) for the macro- 
scopic wave function \T/ (r, t) (also called the order parameter): 



ih 



dt 



2m 



V 2 + V ext + g\V\ 



(1) 



l 



The external potential for a harmonic trap is V ex t(r) = m J2i u i x i /2 and g = 4Trh 2 a s /m characterizes the nonlinearity 
which depends on the particle interaction strength through the scattering length a s . The ground state ^f g is the lowest 
energy eigenstate of the condensate and a solution to the time-independent NLSE: 



2m 



(2) 



where the energy of the ground state fi the chemical potential of the system. 

One way to derive the collective excitations is to linearize the GPE for small perturbations around the ground state 
with the ansatz 



*(r,t) 



-ifit 



(3) 



Substitution into the GPE and linearisation with respect to the small ampitudes hi yields the Bogoliubov-de Gennes 
(BdG) equations 



Ci 



■ g 

,*2 



The operator C is given by 



C 



j-2 

-—S7 2 + V ext (r)-fi + 2g\* g \ 



(4) 



(5) 



By solving the BdG equations we find the eigenmodes with energies hwi, and wavefunctions Vi,Ui that satisfy the 
orthogonality and symmetry relations: 



d 3 r 



(Tr i 



UiU. 



= Si 



= 0. 



(6) 



The small complex amplitude coefficients bi , b* in Eq. (|3|) can be replaced by annihilation and creation operators 
i>i,b\ respectively. This is justified by the standard approach of second quantization, where the eigenmodes of a 
classical system are found and then the complex amplitudes are replaced by mode operators. Alternatively one can 
start with the Grand Canonical many-body Hamiltonian for the field operator ty(r,t), 



H= d 3 r¥(r,t) 



■^V 2 + V ext (r)-[x+^(r,t)*(r,t) 



*(r,i), 



and make the ansatz: 



*(r,t) = tt s + J2 U(r)h(t) + v*(r)b\(t) 



(7) 



(8) 



In this approach the field operator is split into its expectation value (the condensate part) and a fluctuating part 
that accounts for collective excitations and the thermal cloud. Substitution of Eq. (|8j) into the Hamiltonian of Eq.(^|), 
and neglecting terms of order three or four in the excitation operators hi , b\ gives a quadratic Hamiltonian which is 
diagonalized exactly if tp g satisfies the GPE of Eq.(||) and the wavefunctions u.^Vi are solutions of the BdG Eqs.(g). 
The Hamiltonian can therefore be written as 



H = E g + J2 f^Mk + C. 



(9) 



Here C is the zero-point energy of the non-condensate and E g is the energy of the condensate given by 



E„ = / d A r^ 



2m 



(10) 



So far, in all experiments on collective excitations the eigenmodes have been excited strongly into a coherent state. 
For these conditions one can assume that the mode operators commute and replace them by complex numbers so that 
Eq.(||) reduces to Eq.(||) except for a factor e~ v * which amounts to a shift in the zero of energy. 



2 



III. CALCULATING THE QUASPARTICLE WAVEFUNCTIONS IN THE HYDRO-DYNAMIC LIMIT 



In this section we give a brief overview of how to calculate the excited state wavefunctions ttj, V{ directly following 
the approach given in ||. The starting point are the BdG Eqs.(^). These can be rewritten in dimensionless units by 
introducing the following coordinate transforms: yj — — 1..3, where lj = (2fj,/mLu 2 ) 1 ^ 2 are the characteristic 

lengths of the condensate in the Thomas Fermi regime. As in || we define the small dimensionless parameter 
£ = frbj/2/j, and the dimensionless energy of mode i, e$ = Ei/ttio, where Ei is the energy of mode i and ui — {uj\uj2UJ^) 1 ^ 3 ' 
is the geometric mean oscillator frequency. We also introduce the mean characteristic length l c = (hhh) 1 ^ 3 and the 
dimensionless Laplace operator A = ^j = i{wj/to)d 2 /dy 2 and define y 2 = Ylj=i Vj- The resulting equations are: 

- £ 2 Au 4 + y 2 Ul + {2 Ul + v t )n Q = (1 + 2^)^, (11) 
-£ 2 Av, + y 2 v l + (2vi + Ui)n = (1 - 2^)^, (12) 
-£ 2 A?/> 3 + y 2 ip g + n a ip g = i/) g) (13) 

where uq = \ip g \ 2 g/ /i. These equations can be combined to form fourth order equations for the functions ff 1 = Ui±Vi. 
In the hydrodynamic limit the expression ip g — y^no(l — y 2 ), where no = /i/ 1 g is the maximum condensate density, 
is then substituted for the ground state wave function tp gi and terms of second order in the small parameter £ are 
omitted. 

Now we introduce the operator G with the definition: 

G= (l-y 2 )A-2^(^/w) 2 d/%. (14) 

i 

and define new functions Wi(yi, y 2 , 2/3) by ff(y) = Cf(l — y 2 ) Tl ^ 2 Wi(yi,y2,y3), where the relation between the 
coefficients Cf~ is given by Cf — e£C 4 ~. Using these definitions one finally obtains from Eqs.(]ll],|l2|) the compact 
expression Q: 

GW + 2e 2 IU = 0, (15) 

where we have omitted the mode-index i for simplicity. For a spherical trap this equation can be solved exactly. It 
is a hypergeometric differential equation with Jacobi polynomials as the general solution. The quantization of the 
energies comes from the condition that the function must converge at the condensate boundaries which yields an 
analytic expression for the mode spectrum [^|^] . 

In the most general case of an anisotropic trap with three different trap frequencies, one can make a polynomial 
ansatz. The symmetry of the Hamiltonian means that parity is a good quantum number for any spatial coordinate. 
We shall find the solutions for the quadrupolar modes where the order of the polynomial is 2. We can make the ansatz 
§ 

W oc yiy-j, i ^ j (16) 

to find the three odd parity eigenfunctions with eigenfrequencies il = (uof + u 2 ) 1 / 2 . These so-called scissors modes 
have been studied extensively by our group Q and we use Eq.([l6|) to derive some of their coupling properties in 
section [v| We have to make a different ansatz to find the three even parity eigenfunctions (which are also referred to 
as diagonal quadrupolar modes 0): 

3 

W^l + ^^/^fy 2 . (17) 

3=1 

The polynomial coefficients bj completely characterize the mode geometry and will be very important later on in our 
expression for the coupling matrix element. In the following we use the abbreviation bj = bj(u)/ujj) 2 . Substitution of 
Eq.(|l7|) into Eq.([l|) shows that W is a solution provided that the following equations hold 

(18) 
(19) 
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where the matrix S is defined as 



S = 



1 

V 1 



1 

1 



1 
1 

3-^ 



(20) 



The eigcnfrequencies are found by demanding that det(S) = 0. The resulting equation is 
n 6 - 3ft 4 (o; 2 + uj\ + wl) + 8ft 2 (w 2 w 2 + u\ul + ujjujj) - 20(w 2 w§ct|) = 0. 



(21) 



This general expression simplifies for the case of an axially symmetric trap (tui =102). In this case the solutions are 
51 = \f2u]\ for the m = 2 mode, and 



Q 2 = 2 + ^A 2 T^^ 



16A 2 + 16 



(22) 



for the m = low-lying and the m = high- lying mode. Here, the trap anisotropy is given by A = u> z /u> r , where 
uj z ,uj r are the axial and radial trap frequencies respectively. Eqs. ( |2l| , |22| ) were derived in the review paper on BEC 
by Dalfovo et al. In the next step the polynomial coefficients bj are found from any two of the three equations in 
(|f) and Eq.@. 

So far we have summarised important known results that enable us to calculate the quasi-particle wavefunctions for 
all six quadrupole modes. We will use these in section ^ to calculate the matrix elements which describe the couling 
of these modes. 



IV. A MODEL FOR SECOND HARMONIC COUPLING BETWEEN TWO MODES 

We follow the approach given in [Q to derive a set of coupled nonlinear equations (describing second harmonic 
generation) from the NLSE. For convenience we normalize the condensate wavefunction to unity and change the 
parameter g in Eq.([l]) to Nog, where Nq is the number of particles in the condensate. We introduce a set of excitations 
that is normal to the condensate and also diagonalises the many-body Hamiltonian of Eq. (0). This is achieved 
by projecting out the overlap with the condensate from the solutions to the BdG equations to give quasi-particle 
wavefunctions defined by 

Ui = Ui - Ci^g (23) 
v*=v* + c** g , (24) 

where q = J d 3 r [^ g Ui] = —J d 3 r [9 g Vi]. These wavefunctions still diagonalise the many-body Hamiltonian (^) and 
the orthogonality relations (ra) hold as well. The advantage of introducing excitations orthogonal to the ground state 
is that it makes it easier to extract the amplitudes of various excitations from a given wavefunction. In terms of the 
orthogonal excitations a general wavefunction can be written as jij : 



I i>0 ) 



(25) 



where the coefficient b g describes the change in the condensate. It is easy to show that for the orthogonal excitations 
the following relationships hold 

J d 3 r V;*e +V * = l + b g 

J d 3 r [u*^e+^ t - t^e-*"*] = (26) 

The population of the condensate ground state is given by |1 + b g \ 2 No and the population of the excited states by 
\h\ 2 N . 

In the next step we obtain the equations of evolution for the complex coefficients bi(t) by substituting the expansion 
of the wavefunction (25) into the GPE dll), and carrying out the projections described by Eqs. (129) . We so obtain the 
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Heisenberg equations for the c- number equivalents of the mode operators hi. We then transform these equations for 
the mode amplitudes into the interaction picture by making the ansatz = bf-(t)e~ lUit . This gives rise to a large 
number of terms oscillating at frequencies u>i ± u>k ± . If we focus on second harmonic processes where cuk = and 
LOj ~ 2ui we can neglect all the rapidly oscillating terms and retain only the term oscillating at Ay = luj — 2u>i. This 
is called the Rotating Wave Approximation (RWA). If we neglect any variation in the population of the condensate 
mode we obtain the following coupled equations of motion for the two modes i = 1 and j = 2: 

lh ^L = N gM 12 b?*bK e - lA ^ (27) 
at 

ih^j- = \N gM* l2 b^e^\ (28) 

where the matrix element Mi 2 is given by 

Mu = 2 J d 3 r (2ulv*u 2 + iJJvJtia) + % {2u\v{v 2 + u^u 2 )] . (29) 

These equations describe the transfer of excitation between the two modes via annihilation (creation) of two quanta 
in mode 1 and creation (annihilation) of one quantum in mode 2 which is also called a second harmonic process. The 
matrix element M\ 2 contains all the information on the geometry of the two modes that are coupled. If we excite 
the lower mode at resonance (A12 = 0) and there is no initial population in the upper mode, then all the excitation 
is transferred to the upper mode. The opposite is not true, i.e. there is no transfer from an initially excited upper 
mode to the lower mode if we start off with zero population in the lower mode. We will see later in this section 
that in a quantum mechanical description, where the lower mode is described by operators rather than c-numbers, 
downconversion does occur. The strength of the processes depends on the spatial overlap between the respective 
quasiparticle wavefunctions. The characteristic time scale for the transfer from mode 1 to mode 2 is given by the 
expression: 



T = 



V2h 



NgMMO) 



(30) 



This model for the second harmonic coupling between two collective excitations allows us to find an explicit expression 
for the matrix element governing the process. 



A. A Quantum Mechanical Model for the Coupling 

Alternatively, the two coupled nonlinear equations ( p8| ) can be derived from the Hamiltonian (|3l|), which gives a 
full quantum mechanical description and clearly shows the underlying physical processes 

H = fhuia\ai + 2friL>ia\a 2 + — (a\ 2 a 2 + a\a\), (31) 

where a\ai, a 2 a 2 give the quasiparticle populations of mode 1 and 2 respectively. We derive the Heisenberg equations 
for the mode operators by the relation 

&i = % - [H, a % ] . (32) 



To remove the fast oscillation at the mode frequencies lui,u> 2 from the operators a\, a 2 we introduce the slowly varying 
operators b\,b 2 : b\ = ai^==e Ml *, b 2 = a2y=e IW2 *, where the v^o factor arises because of the different normalization 
The Heisenberg equations in terms of these operators are 



ih lT = V^K&i^e-^ 12 ' (33) 



th lT = ^V^o^Mie jAl2t . (34) 



Eqs.(Eq) are obtained by replacing the mode operators by complex numbers and setting 



5 



N Hk = N gM 12 . 



(35) 



However, quantum effects such as squeezing and sub-Poissonian statistics in the quasi-particle number can only be 
described by the operator equations ( |34| ) and are lost in making the classical approximation leading to Eqs.(|2^), where 
operators are replaced by complex numbers. 



B. Squeezing in Parametric Down- and Up-conversion 



In this subsection we want to show briefly how a quantum description of the nonlinear processes leads to nonclassical 
effects such as squeezing. We apply the well established theory of nonlinear effects in optics |}],[l0| to describe the 
phononic coupling in a condensate and demonstrate the direct dependence of squeezing on the coupling matrix 
element in both, up- and down-conversion processes. It is difficult to study squeezing for the full quantum mechanical 
description of both modes given by the Hamiltonian ( j3lj) as the operator equations are nonlinear. We first investigate 
down-conversion for resonant coupling and transform the Hamiltonian ( ^l|) into the interaction picture to obtain 



(36) 



where Si, b\ denote the operators a%, a\ transformed into the interaction picture. Then the mode operators &2, b\ for 
mode 2 are replaced by the c- number (3 2 (we can without loss of generality assume that f3 2 is real), but we retain the 
operators for mode 1. In addition we will treat the mode amplitude of the upper level as a constant. This assumption 
does not account for the depletion of the upper mode and is only valid for small times when the occupation of the 
lower level is much smaller than the occupation of the upper level, i.e. for N\ <C (3 2 . The resulting Hamiltonian is 
quadratic in b\, b\ and gives linear Heisenberg equations: 



dt 



= [h,H R ] = K p 2 b\ 



=£ = [blH R ] = Kfobi. 



(37) 



Eqs.(p7|) can be diagonalised by expressing them in terms of the two quadrature phase amplitudes Q X ,Q P defined 



as 



Qx = h + b 
A h-b 



(38) 



Simple integration yields Q x (t) = e Kl32t Q x (0), Q p (t) — e Kl32t Q p (0). These solutions can be used to calculate the 
evolution of b\ , b\ and the evolution for the number of down-converted quasi-particles, N\ , for which we find (assuming 
mode 2 was initially in a vacuum state |0)): 



Ni = (0|&l6i|0) = smh z {np 2 t). 



(39) 



Eq.(p9[) shows that in this quantum model downconversion occurs even for zero initial population in mode 1. This is 
in contrast to the result of the semiclassical model discussed above. The evolution of the variances in the quadrature 
operators is found to be: 



AQ x (t) 



_ e 2re/3 2 t 



AQp(i) = 



_ -2K(3 2 t 



AQ x (0) 
AQ p (0) 



(40) 



This clearly demonstrates the squeezing in the Q p quadrature component. However, it is important to keep in mind 
that Eqs.(|o|) are only valid for short times before the assumption that the upper mode is not depleted breaks down. A 
possible way to avoid depletion of the upper mode is to keep exciting it. The standard way to do this is to mechanically 
force the condensate into oscillation at the frequency and geometry corresponding to the upper mode. 
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Similarly, squeezing occurs during SHG where two phonons from the lower mode are converted into one phonon of 
the upper mode. In this case we cannot replace the lower mode by a c-number as we have done for our investigation 
of downconversion. We can again try to find an approximate solution valid only for small times. This is demonstrated 
in pi , where a Taylor series expansion is used to describe the time evolution of the mode operators. We assume that 
mode 1 is initially in a coherent state defined by = /3i|/?i) with /?i = l^ile 1 !* and mode 2 is in the vacuum 

state. The result for the squeezing of the quadrature Q x (of mode 1) to second order in time is then given by: 

= 1 - i K 2 t 2 |/3 1 | 2 cos(2$) + 0{gtf, (41) 

This result only holds for small times for which -||/3i| 2 K 2 t 2 <C 1. However, we can see how in both cases, down- and 
up-conversion, the squeezing of the quadrature components can be directly related to the nonlinear coupling strength. 



V. CALCULATING THE NONLINEAR COUPLING RATE 



We will now calculate the nonlinear coupling rates between two condensate excitations and compare them to those 
from recent experiments. For convenience we take linear combinations of the normalized quasiparticle wavefunctions 
to give the new functions , f~ defined by: 

f+ = u t +h = f+ (42) 
fi = Ui - Vi = fr - 2citp g . (43) 

Written in terms of these functions the matrix element in Eq.(p9|) has the form 

Mia = 2 / d 3 r ^ \ ~/+* (/+*/+ + frfc) + ]f+ (/+*/+* /f */f *) } , (44) 



where we assumed that ip g is real. Alternatively, M12 can be written in terms of untilded functions (/j + , f i ) as the 
sum of two parts, M\% = Ad[^ + M^ 2 ) , which are defined as follows: 



M# = 2 J d 3 r ^ g I X -f+* (/+* + f+ + /f*/-) + 1/+ (/+*/+* - f-*fr) } , (45) 
M$ = 2 J d?v ^ g {/+* (2c lC2 ^ 2 - ctV s / 2 - - ca^/r*) + H OWf* ~ <?tf) } • (46) 

(2) 

.M 12 is zero if neither of the quasiparticle wavefunctions has any overlap with the condensate ground state. We can 
now use the functions f^~, /j~ which we found in section III. We can write in general: 

/+ = A iei £(l - y 2 y x ' 2 W x , /f = A 1 {1 - y 2 fl 2 W^ . 
/+ = A 2 e^(l - y 2 )-V 2 W 2 , = A 2 (l - y 2 )^ 2 W 2 , [i '> 

where the A\,A 2 are normalization constants determined from the normalization condition (^[) and 
Wi (yi,y 2 , 2/3)) W 2 (yi, j/2) Va) are solutions to Eq.(|l5|) for mode 1 and 2 respectively. Substituting these expressions 
into Eq.(p5|) gives: 

M$ = J d 3 yW 2 W 2 [3e 2 e 2 £ 2 (l - y 2 )' 1 ~ Ae 12 (l - y 2 )} . (48) 

The first term in the above integral proportional to (1 — y 2 )^ 1 diverges at the condensate boundary (y = 1), but it 
must be dropped as it scales proportional to £ 2 and in the derivation of the quasiparticle wavefunctions we omitted 
terms proportional to £ 2 in the governing equations (hydrodynamic approximation) to obtain Eq.(|l5|) for W(yi, y 2 , ya)- 
We will see later in this paper that this is fully justified by comparison with exact numerical calculations. Note that 
the second term equals zero if the detuning Aei 2 = e 2 — 2t\ is zero. 
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A. Coupling between two even parity quadrupolar excitations 



So far we have made no assumptions about the geometry of the two modes that are coupled and Eqs.(]45j-[16|) are valid 
for any pair of modes. We now focus on the quadrupolar modes of a triaxial trap and investigate the coupling between 
any two diagonal modes, for which the function W{y\, y2, 3/3) is represented by the polynomial given in Eq. (|l7|) . We 
can calculate the matrix element M12 from Eqs.( |46|j48| ). An explicit expression for Mi 2 and its derivation is given 
in Appendix A. It is important to note that for on-resonant coupling between the two modes (u>2 — 2u>i) the matrix 
element simplifies considerably. In that case we can give a simple expression for the radial integrand of the matrix 
element in terms of the dimensionless position y, where the condensate boundaries are given by y = 1: 

""-WI*f ^-^--VsNs- (49) 



where A2 denotes the normalization constant for the wavefunction of mode 2, given by Eq.(AS). It is important to 
note that Eq.(^) describes the resonant coupling between any two diagonal quadrupole modes in a triaxial trap. This 
allows a quick and easy calculation of the coupling strength, coupling rates and squeezing effects associated with the 
nonlinear process. We can see from the radial integrand of Eq.(f49|) that the coupling is strongest in the outer regions 
of the condensate and reaches a maximum at a distance of 0.87 l c from the center. In this region the condensate is 
still well described by the hydrodynamic approximation which only breaks down at a distance of order the healing 
length from the condensate boundaries. The healing length in our recent experiments WM] was about 0.05 l c . 



B. Coupling between two modes in a spherical trap 



Now we want to give a quantitative comparison between the solutions we found to the coupling matrix element in 
a hydrodynamic approach and the exact solutions calculated from the numerical solutions to the BdG equations. To 
facilitate the numerical calculations we will look at the coupling between two quadrupolar modes in a spherical trap 
of frequency to, where the total angular momentum I and the azimuthal angular momentum m are good quantum 
numbers. The two modes are the I = 2, m = mode and the I — 0, m = mode with frequencies of \/2lj and y/5u> 
respectively. In a trap with only axial symmetry these two modes get mixed and become the m = low-lying and 
the m = high-lying (breathing) mode. Their quasiparticle wavefunctions can be presented by the ansatz ( |47j ) where 
W\ — (3cos 2 # — l)y 2 and W2 = (1 — 5y 2 /3). These are the solutions to the hypergeometric differential equations 
discussed in section III . The normalization constants are At = v/35/167rZ3ei£ an d A 2 = (3/2)^/7 /An l^. 

One can see from Eq.(^6|) that in all terms containing the constant c\ are zero, because the overlap between 

mode 1, which is proportional to Y§ , and the ground state, which is proportional to Y °, is zero. The only remaining 
term is proportional to C2 and it turns out to be — 2y^no / NoA2y 6 (1 — y 2 ), which is exactly the expression we found 
for the resonant case (A12 = 0) in a general triaxial trap (see Eq.([49|)). But for these two modes A12 7^ and we 
have to consider the contribution from as well so that we obtain 



M12 = 




— Ai2(l-y 2 )(l-5/3y 2 ) + (l- 2/ 2 ) 
4ei 



y 6 dy 



(50) 



The wave functions for the ground state and the coupling matrix elements M12 are plotted in Fig. ( |T|) a nd Fig.(||) 
respectively. It is important to note that for the resonant case the only contribution comes from Eq.(49). Also, for 
the not quite resonant case displayed in Fig.(||) the integrand is dominated by this contribution. This shows that 
the coupling between different quasi-particle excitations predominantly takes place in the boundary region of the 
condensate and an explicit analytic expression for the spatial probability of the nonlinear process is given by Eq. (E9ft . 
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FIG. 1. The wavefunctions for the ground state in the hydrodynamic limit (solid line) and the exact solution to the GPE 
(dotted line) plotted in units of against the distance from the center of the trap in units of the characteristic length l c . 

The trap is spherical with 1.5 x 10 4 atoms and a frequency lo — 120 Hz . The healing length f; = 0.05Z C 
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FIG. 2. The radial integrands for the matrix-elements M12 in the hydrodynamic limit (solid line) and for the exact numerical 
calculation (dotted line) plotted in units of l~ 3 against the distance from the condensate center in units of l c - The coupling is 
between the I = 0, m = and the I = 2, m = mode for the same trap conditions as in Fig.(|l|). 

The integrated values of Myi in the hydrodynamic approximation and for the exact numerical calculation are 
—0.161 l~ 3 and —0.157 Z~ 3 respectively. The error of the approximate analytical result with respect to the exact 
numerical calculation is in the order of £ 2 , as we expect. The good agreement between the two, even for a relatively 
small number of atoms, justifies the hydrodynamic approximations made in calculating the quasi particle wavefunctions 
and the coupling matrix element. 

C. Coupling between two even parity modes in a TOP-trap 

In our experiments we use a TOP-trap which is axially symmetric and has an anisotropy denned by the parameter 
A = uj z /uj ri where uj z and u) r are the axial and radial trap frequencies respectively. In a recent experiment we studied 
the coupling between the m = low-lying and the m = high-lying mode (which arise from the I = 0, m = and 
the I = 2, m = mode of the spherical trap when A ^ 1) . We can use the formula for M\i for the general triaxial 

) in the Appendix, to calculate the matrix element for any trap geometry. The result 
is shown in Fig.(|3|) in terms of the dimensionless quantity 77112 = MyiI\sJ\. The dependence of A/12 on number and 
mean- frequency is contained in ?,T 3 £ -1 ^ 2 and thus mi2 only depends on the mode geometry, i.e. it is only a function 
of A. In order to get M12 for a specific trap, one has to read the dimensionless matrix element m\% from Fig.(||) and 
multiply it by Zj 3 £~ 1/2 . 





trap, given by Eqs.flAlC . All 
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FIG. 3. The quantity ran — M\il\\/\ is plotted against the trap anisotropy A. The two vertical lines show where there is 
resonant nonlinear coupling between the m = low-lying and the m — high-lying modes. The resonances are determined by 
the matching of mode frequencies (2u)i = UI2) and they are located at A = 0.68 and A = 1.95. 




From Eq.(|22|) one can find two values for the anisotropy A where the two m = modes are resonantly coupled. 
These resonances are at A = 0.68, 1.95 with |m 12 | = 0.038, 0.028 respectively. The characteristic coupling time is given 
by Eq.([50|) which can be written as follows: 

We can calculate the coupling rate for the parameters of our second harmonic generation experiment [Q , which were: 
A = 1.95, Nq = 1.5 x 10 4 , uj r ~ 120 Hz. Inserting the relevant quantities into Eq.(|5l|) we obtain a transfer time of 5.7 
ms and 3.7 ms for small initial excited populations of |&i(0)| 2 = 0.02 and |&i(0)p = 0.05 respectively. For stronger 
initial excitation the transfer times would be even smaller. One can see from Eq.(|5l|) that for larger atom number and 
stiffer traps the coupling times become smaller as well. This shows that the process of coupling the two quadrupole 
modes happens on a timescale of a few milliseconds. It is consistent with our experimental observation Q| where the 
second- harmonic was observed as soon as the driving of the fundamental finished (excitation time of about 30 ms). 



D. Coupling between a scissors mode and another quadrupolar excitation 



Now we will discuss coupling between the off-diagonal quadrupolar excitations. We already showed in section III 



that these modes have odd parity and are characterized by W oc yiUj, i ^ j- Thus their overlap with the condensate 
ground state is zero and the only contribution to the coupling matrix element comes from given in Eq. (^8|) . But 
this integral is zero because the product Wf W2 is odd. So the total matrix element M12 equals zero and there is no 
direct coupling via a second harmonic process between two scissors modes. Our recent experimental observation of 
a downconversion process between two scissors modes must therefore have a more sophisticated interpretation than 
the conversion of one quantum of the higher mode into two quanta of the lower scissors mode. 

Similarly, there is no coupling between a higher-lying scissors mode and a lower-lying diagonal quadrupole mode. 
In this case M$ of Eq.© is zero beacuse W2 is odd. of Eq.(|46]) is also zero because c-i — and f% are odd. 

Thus the total matrix element is zero. However, there is second harmonic coupling from a lower scissors mode to a 
higher lying diagonal quadrupolar mode. For resonant coupling two quanta of the scissors mode are converted into one 
quantum of the higher lying even parity mode and one finds that the matrix element is again given by expression (E^). 



VI. CONCLUSION 



Starting from the NLSE we have derived a simple model describing the nonlinear coupling between two modes. We 
have demonstrated a way how to calculate analytically the matrix element governing the coupling process. We then 
focused on the quadrupole excitations and found that all resonant (u>2 — 2lui) direct coupling processes between the 
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six quadrupole modes are described by the expression in Eq.(^9|) unless they are forbidden (Mi 2 = 0). All second 
harmonic processes involving odd parity (scissors) modes are forbidden except for upconversion from a scissors mode 
to a higher lying even parity mode. It is possible to show that there are other allowed nonlinear processes involving 
all three of the scissors modes. This gives rise to nondegenerate parametric amplification and multimode squeezing. 
Full details and the derivation of an effective Hamiltonian describing all allowed nonlinear processes between the 
quadrupole modes (not just the second harmonic generation described here) will be given in a future publication. 
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APPENDIX A: MATRIX ELEMENTS FOR THE DIAGONAL QUADRUPOLAR MODES 

We want to calculate the coupling matrix element for any two diagonal quadrupole modes in a triaxial trap for 
which the function W(yi, 1/2, J/3) is represented by a polynomial as given in Eq.Jl7|). It is useful to derive a number of 
relations for the polynomial coefficients which allow us to simplify the expressions for normalization constants, overlap 
coefficients and the coupling matrix element. 

The lowest mode of Eq. (||) is the so called Goldstone mode with cjq = 0, tio(r) = v E'g(r) and Vo(r) = — ^*(r), 
which arises from the U(l) symmetry breaking. For this particular mode the orthogonality and symmetry relations 
of Eqs.(^) take the form 

J d 3 r(V* g u t + y g v t ) -0. (Al) 



Eq.(Al) implies that J d 3 Ttp g f + — 0. Substituting (|47| ) into this equation and integrating over the angles and radial 
coordinate gives: 

£> = -5 (A2) 

3 

The characteristic polynomials for the quadrupole modes are real and thus Uj , V\ can be taken as real which allows us 
to derive from the orthogonality relations (g): 

' >x/+/r (A3) 

If we now insert two different polynomials corresponding to two different solutions for f^, jj into Eq.( |A3|) we obtain 
the relation J cPyWiWs = and from that: 

X]Mi= 5 (A4) 

i 

M2 + Mi + hd 3 + hdi + M3 + M2 = 20, (A5) 

where the fej's denote the polynomial coefficients of mode 1 and the d^s the polynomial coefficients of mode 2. These 
coefficients are found from Eqs.([l8|,[l9|) in section H. We also have to calculate the overlap coefficients Cj between the 
condensate and the quasiparticle wavefunctions: 

d 3 r^>, = \[ d 3 r^I(/+ +fr) = l [ d 3 rfr. (A6) 



2 / TgyJl Jl ' 2 



After substituting Eqs.(|47|) into Eq.(A6) we obtain 



1 . ]3 /n^87T 

2 1 C V N Q 15 



(l^E^=^, (A7) 
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where we used (A2) and No = nol^8Tr/15 (Thomas Fermi relation for fi(No)) for the last step. We also have to 
calculate the normalization amplitudes for these modes and obtain from Eqs.([A3|): 



1 =l««iffi 



3 ~ h ) + 2 (^i^2 + Ma + hh) ~ 35 



-1/2 



(A8) 



We can now calculate the coupling matrix element from Eqs.(|4£],[4^). We shall introduce a constant R\ 2 defined as 
follows: 



J?12 



The first part of the matrix element is then: 



M$ = -R12 J dy{ ,/ 



no -A\A^ll/2 



47T 

105' 



(A9) 



15 b ldi + 3(2&iMi + b\d 2 + 2biMi + bjd 3 + 2bib 2 d 2 + b\d Y 

i 

26i6 3 J 3 + &pi + 26 2 6 3 J 2 + b\d 3 + 2b 2 b 3 d 3 + b 2 3 d 2 ) + 2b 1 b 2 d 3 + 2bib 3 d 2 + 26 2 6 3 di 
490 + 21 ^2 b\ + 7(26x62 + 26 2 6 3 + 26361) 



+ 



- 525^ + 105^^12(1-2/^2/. 



(A10) 



The second part of the matrix element is given in Eq.(All). Note that this part arises due to the finite overlap 
between the condensate ground state and the untilded Bogoliubov wavefunctions Ui,Vi. 



A4 2) = AR 12 / Ae 12 I 352/ 4 - + ^ \ y l {\ - y*)dy 



7 



-2J^A 2 / y\l-y 2 ).dy 
iV o Jo 
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(AH) 



We note that for resonant interaction Aei 2 = e 2 — 2ei = 0, so that = and the only term surviving in M^ 2 ) is 



-2^ ~/N~ A 2 J Q 1 y 6 (l~y 2 )dy. 
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